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Path integral Monte Carlo (PIMC) simulations are used to calculate the momentum distribution 
of the homogeneous electron gas at finite temperature. This is done by calculating the oft-diagonal 
elements of the real-space density matrix, represented in PIMC by open paths. It is demonstrated 
how the restricted path integral Monte Carlo methods can be extended in order to deal with open 
paths in fermionic systems where a sign problem is present. The computed momentum distribu- 
tion shows significant deviations for strong correlation from free fermion results but agrees with 
predictions from variational methods. 



I. INTRODUCTION 

The momentum distribution is one of the fundamental 
properties of a quantum system. It can be directly mea- 
sured by inelastic scattering or by studying the trajecto- 
ries of particles. The first method has been used to deter- 
mine the condensate fraction in liquid helium while the 
latter was used to demonstrate that Bose-Einstein con- 
densation of supercooled alkali atoms in magnetic traps 
had been achieved. 

In this article, we calculate the momentum distribu- 
tion of fermion systems at finite temperature with path 
integral Monte Carlo (PIMC). The results show how the 
fermionic momentum distribution evolves as a function 
of temperature. At a temperature much higher than 
the Fermi temperature, Tp, the fluid has a Maxwell- 
Boltzmann momentum distribution and with decreasing 
temperature, the Pauli exclusion principle becomes more 
important giving rise to a Fermi surface. While it is very 
simple to study this cross-over in noninteracting systems, 
interparticle interactions require a much more sophisti- 
cated description. 

Here we compute the effect of interactions of the mo- 
mentum distribution at finite temperature using fermion 
restricted path integral Monte Carlo simulations. To de- 
termine the needed off-diagonal density matrix elements 
we need to perform simulations with open paths. We 
compute these properties for the homogeneous electron 
gas and demonstrate how increased correlation effects al- 
ter the momentum distribution at different temperatures. 



II. PATH INTEGRAL MONTE CARLO 

We now give a brief review of the PIMC simulation 
technique for fermions and of the restricted path method. 
In PIMC calculations of the diagonal density matrix 
each particle is represented by a closed path in imagi- 
nary time. Following the procedure developed for bosonic 
systems Q, we extend fermionic PIMC to estimate off- 
diagonal density matrix elements needed for the momen- 
tum distribution. This requires one of the particle paths 
to be open. Special emphasis will be placed on how the 
path restriction is applied to open paths and how the 
Monte Carlo sampling is affected. 



A. Restricted paths technique 

The thermodynamic properties of a quantum many- 
body system can be derived from the thermal density, 
p = e~ l3n with (3 = 1/ksT. For the purpose of perform- 
ing Monte Carlo simulations, we express this operator in 
position-space, 



p(R,R';/3) 



(R 



!*'> = £< 



-pe. 



*:(r) * s (r') 



(2.1) 

where ^ s are the many-body eigenfunctions and e s the 
corresponding eigenvalues. R = {ri, . . . , r^} represents 
a set of coordinates of N particles in d dimensions. The 
density matrix of a bosonic (B) or fermionic (F) system, 
/Ob/f(Ri R'; P), can be constructed from the density ma- 
trix for distinguishable particles by a sum of permuta- 
tions, V, to project out states of the corresponding sym- 
metry, 
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where (±l) p denotes the sign of the permutation. Using 
the operator identity, e^^ u = (e~ TH ) M , the density ma- 
trix at temperature T can be expressed in terms of den- 
sity matrices at a higher temperature MT. This leads to 
a path integral with M steps in imaginary time of size 
r = (3/M, 

Pb /f(R, R'; P) = ^ ^ V j- ■ / dRl dR2 dRM 1 



in a few cases, e.g. for noninteracting particles. In prac- 
tice, one introduces a trial density matrix pt(R, R'>/3) 
that provides approximate fermion nodes which intro- 
duces an approximation in computed observables. How- 
ever, for many systems this technique has worked well. 

The simplest approximation for the trial density ma- 
trix pt(R, R'; 0) is a Slater determinant of single particle 
density matrices, 



pd(R,Ri;t) pd(Ri,R2;t) . . . / o d (Rm-i,£ , R';i~) 



fEW / dR * e 



-S[Rt] 



R^PR' 



where 5 represents the "action" of the path R t begin- 
ning at R and ending at PR'. Here we use the pair 
density matrix for the link actions: /9d(Ri, R2; t). For 
bosonic many-body systems, the integrand is nonnega- 
tive and this expression can be efficiently evaluated us- 
ing Monte Carlo techniques |]| . In the case of fermions, a 
straightforward evaluation of this expression is impracti- 
cal because the cancellation of many positive and neg- 
ative terms of the same order leading to numerically 
inefficient computation in the low temperature, many- 
fermion limit. While one can still use this expression 
to numerically study systems of a few fermions, the effi- 
ciency rapidly decays with increasing number of particles 
and decreasing temperature, ~ e _/3Jv . This is referred to 
as the fermion sign problem. 

Ceperley 0, 0, E] has shown that the fermion sign 
problem in imaginary time path integrals can be solved 
by restricting the path integration to a subvolume of 
the entire path space. Let us define the nodes of the 
fermion many-body density matrix as the surface where 
Pf(R, R'; t) = 0. The nodes are used to confine the paths 
R(t) to regions where the density matrix is nonzero, 
Pf(R*, R(t); t) 7^ 0. R* is called the reference point and 
defines the region in {R, t} space (T(R*;i)) where the 
path is allowed to be. The fermionic density matrix is 
then given by the restricted path integral, 

p F (R*,K';[3) = ±Y,(- 1 f J dR * e_S[RtI > 

V R*^PR'eT(R':i) 

where one sums and integrates over all paths that never 
cross the nodal boundaries. By introducing this restric- 
tion, one effectively cancels all negative and some posi- 
tive contributions to the trace of density matrix. Com- 
plete cancellation of all negative terms, however, is only 
reached on the diagonal of the density matrix. Off the 
diagonal, negative contributions also enter the restricted 
path integral, as will be illustrated below. In either case, 
the restriction gives rise to an efficient numerical algo- 
rithm that scales favorably with increasing number of 
particles, similar to that for bosons. 

The expression in Ea. l2.4l is exact as long as the restric- 
tion is exact ■ The exact density matrix is only known 
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For a homogeneous Fermi liquid, such as the electron gas, 
the single particle density matrices are: 



p{ ) 11 (r,r'; / 3) = (4^A/3) 
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(2.6) 



where A = fi 2 /2m. For temperatures above the Fermi 
energy, where exchange is small, the nodal surfaces of 
this determinant are accurate. At decreasing tempera- 
ture, interaction effects become more important and con- 
sequently the error from employing free particle nodes in- 
creases. (Note, however, that the nodes are second order 
in the interaction.) One can improve the nodal approxi- 
mation by using the dual-reference point method In 
this approach, the nodal restriction is pt(R(£), R*; t*) ^ 
where t* = min(i, (3 — t) is the smaller imaginary time 
separation from the reference point. The dual-reference 
point method will be used throughout this work. 

How sensitive the derived results are to the accuracy 
of the nodes can depend on the interactions in a particu- 
lar system. Generally, one expects free particle nodes to 
work well at high temperature and when correlation ef- 
fects are weak. Also, when particles are localized like the 
electrons in molecular hydrogen or particles in a Wigner 
crystal, the effect of the nodes is reduced and nodal re- 
striction is unimportant. One can improve nodal sur- 
faces in several ways. For example, one can derive the 
nodes from a variational density matrix Q, by treating 
the interactions with Hartree-Fock. Or, one can use back- 
flow in the density matrix to account for interaction ef- 
fects 00. 

Consider the calculation of an observable, O, diagonal 
in R space such as the pressure, kinetic, potential and 
internal energy as well as pair correlation functions: 

(O) = | JdRp(K,R;p) (R|0|R) , (2.7) 

Z = y"dRp(R,R;/3) . (2.8) 

Such computations require only simulations with closed 
paths, beginning at a point R and ending at its permu- 
tation, PR. The restriction eliminates all contributions 
from odd permutations that would enter with a nega- 
tive weight because their paths would violate the nodal 
constraint an odd number of times. 
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B. Computation of the momentum distribution 

The single-particle momentum distribution for N a par- 
ticles in spin state a is defined as, 



(2.9) 



with the normalization for a finite system and in the ther- 
modynamic limit respectively given by, 

E n ( k ) = N " . J dk "( k ) = N ° • ( 2 - 10 ) 

Inserting complete sets of states and using (R|P) = 
e- iR - p / R /(27rn) Jvd / 2 , one finds, 



the observed separations of the open ends weighted with 
the sign of the permutation. The discussion in the previ- 
ous subsection concerning the elimination of odd permu- 
tations holds for diagonal paths. Off-diagonal paths with 
odd permutations do contribute to the Monte Carlo aver- 
ages with a negative weight. At separations where n(s) is 
negative, odd permutations outweigh even permutations. 
The algebraic decay of n(s) requires long exchange cy- 
cles, of the order of the number of particles. In restricted 
PIMC, there is a direct relation Q between long exchange 
cycles and the discontinuity of n(k) at k = kp. 



C. Example of two free fermions 



l(k) = 



zn 



tP-(R-R')/n N " 
dRdR'dP (R|p|R) fn _^ (N _ 1)d £ *(Pi 



zn 



1 JdKdr'^'^pir,, 



,r N ,r' 1 ,T 2 , ■ 



• ,rjv) 



The last line follows from the equivalence of all the par- 
ticles assuming particle 1 has spin a. Consequently, n(k) 
is given by the Fourier transform of the single-particle 
reduced density matrix (SPRDM), n(s), 

n(k) = ^ J ds e-* ks n(s) , (2.11) 

n(s) = J dRp(ri,r 2 , . . . ,rjv , ri + s,r 2 , . . . ,r^) . 

Note that n(s = 0) = 1 from the definition. 

Classical particles have a Maxwellian momentum dis- 
tribution with 

n(k) = ^(4^A/3) d/2 exp{-/3Ak 2 } , (2.12) 



n(s) 



exp 



4A/3 J 



(2.13) 



For an ideal Fermi gas in 3 dimensions at T = 0, the 
momentum distribution (for one spin state) is a Fermi 
function, 

s / 1 for fc < fc F with fc F = (67r 2 iV CT /r2) 1 /3 
?H j ~ \ for k > k F 

n(s) = 3/x 3 [sin a; — x cos a;] with x = skp 

Hence, the free fermion SPRDM decays to zero as 
cos(s kp) / s 2 at zero temperature. 

The function n(s) can be computed from PIMC sim- 
ulations with one open path, where the vector s is the 
separation of the two open ends. Methods developed for 
bosonic systems p], such as liquid 4 He, carry over directly 
to fermion path integrals except for the considerations 
having to do with the nodal restrictions. In the simula- 
tions, n(s) is obtained as a histogram, in which one enters 
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FIG. 1: Illustration of the three different types of paths 
for two free fermions all starting at a given reference point 
(r*,r2). A node-avoiding ("NA"), a node-crossing ("NX") 
as well as a permuting and consequently node-crossing path 
("PNX") are shown in the (ri, r2) plane. The thick solid line 
indicates the node. 

In order to illustrate the restricted path integral tech- 
nique, we discuss the simplest fermionic system: the case 
of two free fermions with the same spin without bound- 
aries. The density matrix px(ri, r 2 , r*, r 2 ; /3) fEa. l2.5|) is 
positive when 



(n - r 2 ) • (ri - r*) > 0. 



(2.14) 



The nodal surface is the hyperplane ri = r 2 at all tem- 
peratures as illustrated in Fig. ^ This figure also shows 
the three types of paths that contribute to the trace of 
the density matrix. Since we only consider diagonal den- 
sity matrix elements in this figure, the path starting at 
the reference point, (r 1; r 2 ), can either return to its origin 
or to the only possible permutation of it, (r 2 , r^). 

In restricted PIMC, only node-avoiding paths, labeled 
"NA" , contribute. In the case of two particles, the re- 
striction prohibits any permutation. In case of the direct 
(unrestricted) fermion method, nonpermuting path that 
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FIG. 2: Illustration of node-avoiding paths that contribute to 
off-diagonal density matrix elements for two free fermions. All 
start at the reference point (r*, r^). Shown are paths with one 
open end at ri. The dot-dashed line indicates possible end 
points for nonpermuting paths [n = r*, labeled "NA"]. In 
contrast to on-diagonal matrix elements shown in Fig.0 here, 
we find permuting but node-avoiding paths with end points 
on the dashed line [ri = v"^,, labeled "PNA"] that represent 
negative contributions to the density matrix. 



cross the nodes ("NX") and those that avoid it, enter 
with a positive sign. Permuting paths ("PNX") now also 
enter with a negative weight because of the (— l) v factor. 
For bosons, permuting and nonpermuting paths enter a 
with positive weight. 

In order to compute the momentum distribution for 
two particles, one would perform simulations with one 
open and one closed path. Fig. [21 illustrates such paths 
in the (ri,r2) plane. For simplicity, only node-avoiding 
paths that begin at (r*,!^) and have an open end at Y2 
are shown. Under these restrictions, only two types of 
paths remain: (i) Nonpermuting paths with ri = r\ and 
Y2 7^ ?2- Path 1 would appear as a closed polymer while 
path 2 would be open, (ii) A new category appears: 
node-avoiding but permuting paths with ri = and 
Y2 7^ r*. In the case of closed paths, the nodal restriction 
eliminates permutations because the final point (rjjS , ) 
lies on the other side of the node. However, for open 
paths, many final points for permuting paths are possible 
(see dashed line in Fig. |5J) because path 2 does not have 
to end at the beginning of path 1. 

With increasing number of fermions, it is difficult to 
illustrate the nodal constraint graphically. However, this 
is not a computational difficulty since one can easily ver- 
ify the sign of the Slater determinant for any proposed 
path to see if it is node-avoiding. 



D. Monte Carlo sampling with open paths 

To apply the restricted path integral to open paths, 
one first has to decide at which time slice one opens the 
paths, relative to the reference point, assumed to be at 
t = 0. Using the dual reference point method, one can 
put the open ends at t = 0, or at t = (3/2. We used the 
latter choice because it is more symmetrical and it avoids 
the singular behavior of the fermion density matrices as 

To move the open and closed paths in the presence of 
nodal constraints, one first picks a time interval of size I 1 
in which the paths will be modified where / is the number 
of levels in the bisection method [lj . If it is chosen too 
large, most trial moves will be rejected. If it is too small, 
the paths diffuse too slowly. The optimal choice depends 
on strength of the interactions and on the fermion de- 
generacy. After selecting the time interval, one uses the 
heat-bath algorithm (see [l| for details) to sample a par- 
ticular cyclic permutation of moving particles. Then one 
proceeds with a bisection method for all I levels in order 
to determine the positions of the moving particles at each 
time slice. In this study, the sampling probabilities are 
taken from the free particle density matrix. For a closed 
path connecting points ri and r2, this becomes, 

TKr) = (2W) d / 2 exp|-^i=l!} , (2.15) 

where r m is the mid point (ri + r2)/2. For strongly in- 
teracting systems, a better sampling efficiency has been 
reached by adding a drift and a covariance term . For 
an open end, which is only connected to one point ri, the 
free particle sampling probability is, 

T,(r) = (4 i 7rAr) d / 2 exp|-^^} . (2.16) 

The pair action, u(R, R';/3) = 

— log[p(R, R';r)/po(R., R';t)], is then used to de- 
termine whether a particular configuration will be 
rejected or temporarily accepted. 

As a final step, one verifies if the nodal constraint is 
satisfied at each time slice. Otherwise, the proposed con- 
figuration is rejected. For fermion simulations with only 
closed paths, one can eliminate odd permutation moves 
because they would inevitably violate the nodes. In the 
presence of open paths, this is not true; otherwise we 
would not get the negative pieces of the density matrix. 
The node-avoiding condition for each time slice is still the 
same, pT(R(t), R*; min(i, — t)) > 0. The following two 
conditions can be used to eliminate paths violating the 
nodal constraint: (a) It is impossible to permute an even 
number of closed paths while keeping all other particle 
coordinates fixed, (b) It is impossible to permute an open 
and a closed path in a move that does not change the time 
slice with the open ends. It should be noted that the rules 
are only employed to improve the efficiency. Paths vio- 
lating these conditions would be rejected anyway when 
the nodal constraints are enforced. 
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TABLE I: Kinetic and potential energy, K and V, per particle 
in units of Hartrees are tabulated for r s — 4 and 40 for the spin 
polarized (JVf = 0) and unpolarized homogeneous electron gas 

(iV T =iV i ). 



2e A t iVj. 


K[r a = 4] 


V[r. = 4] 


K[r s = 40] 


V[r s = 40] 


1 33 

2 33 
4 33 
8 33 
16 33 


0.2961(6) 

0.177(1) 

0.144(4) 

0.143(8) 

0.14(1) 


-0.1513(1) 
-0.16413(6) 
-0.16997(5) 
-0.17223(6) 
-0.17341(4) 


0.00354(2) 

0.002605(6) 

0.00237(1) 

0.00232(3) 

0.00227(4) 


-0.019553(4) 
-0.019774(2) 
-0.019826(2) 
-0.019834(2) 
-0.019834(4) 


1 33 33 

2 33 33 
4 33 33 
8 33 33 
16 33 33 


0.1930(5) 

0.123(1) 

0.105(3) 

0.101(6) 

0.100(8) 


-0.15574(9) 

-0.1627(1) 

-0.1636(1) 

-0.1632(2) 

-0.1631(2) 


0.002859(4) 

0.002389(5) 

0.00230(1) 

0.00227(3) 

0.00232(5) 


-0.019715(2) 
-0.019789(3) 
-0.019803(2) 
-0.019803(2) 
-0.019797(6) 



The SPRDM is proportional to distribution of s, the 
separation of the two ends. Configurations with an odd 
permutation contribute negatively. For a finite homoge- 
neous system, n(s) is only a function of |s|. However, in a 
periodic but finite simulation cell, there is a dependence 
on the orientation of s with respect to the cell boundaries. 
In all following results, n(s) will be spherically averaged. 
However, for the computation of the momentum distri- 
bution, n(k) oc (e lks ), the angular dependence of n(s) is 
important. We compute n(k) directly during the Monte 
Carlo simulation for a reasonable number of k vectors, 
to avoid storing the d dimensional function n(s). 

The normalization of n(s), is not determined a pri- 
ori with a single PIMC run. We use the same method 
developed for the bosonic superfluid helium. The nor- 
malization is determined by its value at n(s = 0) = 1. 
A difficulty arises from the fact that the probability that 
the two ends have a distance less than s scales as s d in d 
dimensions and error bars accordingly increase for small 
s as s~ d / 2 . To smooth the errors, we fit the observed 
histogram of end-to-end separations to a function of the 
form: 



lim nuc(s) = £i 

s— >0 



(2.17) 



for the unknown parameters £i and £2 ■ Here K is the ki- 
netic energy derived from a separate simulation with only 
closed paths. The specific values are given in Table [I]to- 
gether with potential energy, V . The internal energy can 
be obtained by adding K + V and pressure for coulomb 
systems, follows from the virial theorem, 3PS1 = 2K + V. 
The fit then determines the normalization constant £1 of 
the momentum distribution. To further reduce the error 
bars of n(s), we enhanced the frequency of moving the 
open ends by selecting the open path more often and also 
by moving the slice with the open ends more frequently. 



III. RESULTS 

We determined the momentum distribution of the ho- 
mogeneous electron gas at three different conditions. 




FIG. 3: The left panels show the momentum distributions 
n{k) from PIMC simulations (circles) with 33 spin polar- 
ized electrons for different temperatures at (r a = 4). For 
T > Tp/4, the population of low momentum states is en- 
hanced compared with the corresponding ideal Fermion re- 
sults (dashed lines), which leads to a lowering of kinetic en- 
ergy (K < Kq). The right panel shows the corresponding 
off-diagonal density matrix elements n[r). With decreasing 
temperature, a negative region near sIcf ~ 5.8 is apparent. 



First, we study the spin polarized electron gas at a den- 
sity of r s = 4, which corresponds to the electron den- 
sity of a low density metal such as sodium. The com- 
puted momentum distribution turns out to be similar to 
an ideal Fermi gas. Secondly, we calculated n(k) for spin 
polarized electron gas at a much lower density of r s = 40, 
where correlation effects are very large and significant de- 
viations from free particle behavior are found. Finally, we 
present results for the unpolarized electron gas at r s — 40 
where correlation effects are even more important. For 
all three conditions, we have computed the momentum 
distribution for a series of different temperatures ranging 
from 1 < Tp/T < 16. We compare with results from zero 
temperature quantum Monte Carlo simulations. 

Fig. G] shows the momentum distribution for spin po- 
larized electron g function of temperature. In 
the limit of high temperature T ^S> Tp, fermion effects 
become less important and one recovers the classical 
Maxwcll-Boltzmann distribution. With decreasing tem- 
perature, one finds that population of low momentum 
states with k < kp increases until it reaches 1, the max- 
imum allowed by the Pauli exclusion principle. Simul- 
taneously, the slope at the Fermi wave vector becomes 
increasingly steep. We do not exactly recover the limit 
of ideal Fermi function in Eq. 12.141 because of finite size 
effects. The system size we used, N = 33 at T = Tp/16, 
is already a demanding computation, in part because we 
used a time step of r = l/32Tp needed to enforce the 
nodal constraint accurately along the paths but requir- 
ing simulations with up to 384 time slices. Notice that at 




2 4 6 8 10 5 10 15 20 25 30 

s k F cycle length / 



FIG. 4: The left column shows the distribution of the open 
ends, n(s) multiplied by the volume element s 2 (solid lines) 
for a system of 33 spin polarized electrons for r a = 4 at differ- 
ent temperatures (see corresponding Fig.EJ. The dash-dotted 
lines indicate the distribution of paths that enter with a neg- 
ative sign. Their contribution vanishes more quickly for small 
s because they do not contribute to diagonal density matrix 
elements. The dashed lines indicate the distribution of posi- 
tive contributions. The right column shows the distribution 
of permutation cycles as a function of cycle length I times the 
weight of the overall permutation. 

T = TV/16 there are still small but nonnegligible thermal 
excitations of states above kp present. The correlation 
effects, absent for free particles, are small but neverthe- 
less significant. At high temperature, the interactions 
lead to an increased population of low momentum states 
resulting in a lowering of the total kinetic energy. The 
reason for this effect is that the entropy is the dominant 
part of the free energy at high temperature. Interactions 
can lower the entropy which also leads to a lowering of the 
kinetic energy. This effect has been discussed in detail 
in 0. At low temperature, the free energy is dominated 
by the interaction term and the kinetic energy is always 
higher than the corresponding ideal value. As a result, 
even at T = 0, states above the kp are populated. Ac- 
cording to Migdal's theorem 0, as long as the system 
remains a Fermi liquid, the discontinuity at /cf remains, 
but the step size is reduced compared to the free fcrmion 
value. 

We have shown the SPRDM n(s) on the right side of 
Fig. [3J they are very close to those of free particles at 
this density. However, in the more correlated systems, 
discussed later, significant deviations from the ideal be- 
havior are found. At high temperature, (in the classical 
limit) n(s) is dominated by a single Gaussian. With de- 
creasing temperature, a shallow negative region develops 
around skp = 5.8. At skp = 7.25, n(s) becomes posi- 
tive again and exhibits a maximum at skp = 8. Further 
oscillations cannot be identified for this system size. 




k / k F s k F 



FIG. 5: Momentum distribution n(k) and off-diagonal den- 
sity matrix n(r) for the spin polarized electron gas at r s — 40 
from simulations with 33 particles. Deviations from the ideal 
Fermion results are increased compared to Fig. |H] (see descrip- 
tion and line styles there) due to increased correlation effects. 




k / k F s k F 

FIG. 6: Momentum distribution n(k) and off-diagonal density 
matrix n(r) function for the unpolarized electron gas at r a = 
40 from simulations with 66 particles. Correlation effects are 
enhanced compared to Figs. |3]and|S] (see description and line 
styles there). 



Fig. 0] shows the distribution of positive and negative 
contribution of the SPRDM. (The contribution to s 2 n(s) 
is shown.) At small separation, n(s) is dominated by 
the positive contribution because negative terms are pre- 
vented by the restriction since the paths become diagonal 
in this limit. The negative region near skp = 5.8 develops 
because two particle permutations occur with increasing 
probability for temperatures T < Tp. Such open paths 
can spread out further than single open paths and there- 
fore start to dominate at larger s. As the temperature 
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FIG. 7: The temperature dependence of the population of the 
zero momentum state, n(k = 0), is shown for two densities 
(r s — 4 and 40) of the fully spin polarized electron gas (from 
PIMC simulations with 33 particles) and for the nonpolarized 
case (simulations with 66 particles). 



is decreased further below Tp, longer and longer per- 
mutation cycles contribute to n(s). The magnitude of 
positive and negative contributions increases with s but 
each function dominates for different separations giving 
rise to the oscillatory behavior of the total n(s) function, 
which is expected from the zero temperature result. 

Fig. 0] also shows the probability of finding a permuta- 
tion cycle times the overall sign as a function of the cy- 
cle length. At high temperature the thermal De Broglie 
wave length, A^ h = 47rA/3, is short compared to the inter- 
particle spacing, which makes the longer permutation 
cycles occur with small probability. With decreasing 
temperature, longer permutations occur more frequently 
and the cycle distribution approaches a positive constant. 
The occurrence of a particular cycle length is no longer 
correlated with total permutation, and overall there are 
more positive than negative total permutations. The dif- 
ference gets smaller as the negative contributions become 
more important at an even lower temperature. This dis- 
tribution shows oscillation for the largest occurring cycle 
lengths since if such a long cycle occurs it is unlikely there 
is also another permutation cycle that can alter the total 
sign of the permutation. 

Fig. [3] shows the n(k) and n(s) functions for the spin 
polarized electron gas at a much lower density of r s = 40. 
Under these conditions, the particles are significantly 
more correlated and their behavior differs substantially 
from free fermions. As a result, the momentum distribu- 
tion is much more spread out leading to the population 
of higher momentum states. The total kinetic energy 
is clearly above the corresponding ideal value. In the 
low temperature limit, the population of the k — state 
reaches a value of only about 0.8 relative to free parti- 



cles. The n(s) distribution is shifted to smaller s values 
indicating that the paths are more localized due to the 
repulsive interactions. 

Fig. shows the corresponding results for the unpo- 
larized electron gas with 66 electrons. These simulations 
are comparatively easier than simulations at higher den- 
sity because the interactions lead to some localization 
and cut down on the number of long permutation cycles 
which give a fluctuating sign. Simulations of a correlated 
system are therefore less demanding than that of weakly 
interacting particles. 
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FIG. 8: Comparison of the momentum distribution for r a = 
40 at T = Tf/16 computed with PIMC, VDM and ground 
state diffusion Monte Carlo (DMC) calculations at T = 0. 
The upper graph shows results from simulation with 33 spin 
polarized electrons, the lower graph represents the unpolar- 
ized case with 66 particles. 

In Fig. El the occupation of the zero momentum state 
is plotted as a function of temperature. At high tem- 
perature, all curves converge to the free particle result. 
Simulations of the spin polarized electron gas at r s = 4 
are above the ideal result for high temperatures under- 
lining the lowering of the kinetic energy. For low temper- 
atures, they converge to the ideal value of 1 within the 
error bars. The graph also shows that the low density 
results converge to a ground state limit as well. 

In Fig. |SJ we compare PIMC results and variational 
density matrix calculations (VDM) [j| at T = Tp/16 
with the ground state momentum distribution derived 
from diffusion Monte Carlo (DMC) simulations using 
backflow nodes (For additional DMC results see [3). 
The agreement between PIMC and DMC is excellent for 
both spin polarizations considered here. Only for the 
nonpolarized case, one observes some very small devia- 
tions around k ss kp , an indication of thermal excitation 
present. Fig.|§|for the spin polarized case also shows good 
agreement with a VDM calculation using the free parti- 
cle density matrix (Eq. 12. 5|) multiplied by a temperature 
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FIG. 9: Finite size study of the momentum distribution 
for spin polarized electron gas at r s = 40 and Tf/T = 8 
computed using the variational density matrix method for 
N = 33, 57, and 117 particles in periodic boundary condi- 
tions. 

dependent Jastrow factor derived from the random phase 
approximation 8]. However, for the unpolarized system, 
VDM predicts n(k) values that are significantly too high 
for k < kp. While the VDM method can certainly be im- 
proved by choosing parameters in the Jastrow factor Q 
more appropriately, it underlines the need for methods 
that work without input from analytical calculations. 

To estimate the finite size effects, we used the VDM 
method which is significantly less computationally de- 
manding than PIMC. Fig. |5| shows the finite size depen- 
dence of the momentum distribution for the spin polar- 
ized electron gas at r s = 40 and TpjT = 8, a density at 
which PIMC and VDM agree well. The simulations were 
performed for different numbers of particles correspond- 



ing to filled k shell structures with N — 33, 57, and 117 
particles. Using periodic boundary conditions, the func- 
tion n{k) can only be computed for k values in the recip- 
rocal lattice of the simulation cell. Consequently, n(k) 
is shown for a different set of k values depending on the 
number of particles. The overall agreement of the com- 
puted momentum distribution is very good, indicating 
that the finite size errors are small at these temperatures. 



IV. CONCLUSIONS 

This computational technique allows one to calculate 
the momentum distribution within PIMC for fermion 
systems. It combines the sampling methods using open 
paths developed for bosonic liquids with restricted path 
technique derived for fermions. Results for the homo- 
geneous electron gas show that the temperature depen- 
dence of the momentum distribution can be studied and 
ground state results can be reproduced. The method 
is applicable to any Fermi system, in particular to hot 
dense hydrogen |ll| where one expects significant changes 
in the momentum distribution with increasing density 
as the electrons are delocalized in the molecular-metallic 
transition, and to calculate the momentum distribution 
of 3 He [l2( and 3 He- 4 He mixtures. 
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